gusucode.com > 《MATLAB图像与视频处理实用案例详解》代码 > 《MATLAB图像与视频处理实用案例详解》代码/第 19 章 基于语音识别的信号灯图像模拟控制技术/voicebox/momfilt.m
function [y,mm]=momfilt(x,r,w,m) %MOMFILT calculates moments of a signal using a sliding window Y=(X,R,W,M) % % Inputs: x is the input signal % r is a list of moments to calculate % (+ve = relative to mean, -ve = relative to zero) % w is the window (or just the length of a Hamming window) % Note: If the window is asymmetric, you should be aware that it gets % flipped in the convolution process % m is the sample of w to use as the centre [default=ceil(length(w)/2+0.5)] % % mm the actual value of m used. Output point y(i) is based on x(i+m-w:i+m-1). % % Example: % To calculate a running kurtosis using a Hamming window of length 30: % y=momfilt(x,[2 4],30); k=y(:,2)./y(:,1).^2 % Copyright (C) Mike Brookes 2007 % Version: $Id: momfilt.m,v 1.1 2007/07/18 17:42:13 dmb Exp $ % % VOICEBOX is a MATLAB toolbox for speech processing. % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You can obtain a copy of the GNU General Public License from % http://www.gnu.org/copyleft/gpl.html or by writing to % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if nargin < 3 w=hamming(length(x)); elseif length(w)==1 w=hamming(w); end lw=length(w); w=w(:); % force to a column vector if nargin < 4 m=(1+lw)/2; end m=max(round(m),1); mm=m; r=round(r(:))'; % force integer row vector of moments lx=prod(size(x)); lxx=lx+m-1; xx=zeros(lxx,1); xx(1:lx)=x(:); % extend with zeros so filter() works correctly cw=cumsum(w); sw=cw(end); y0=repmat(sw,lxx,1); lxw=min(lxx,lw); y0(1:lxw)=cw(1:lxw); y0(lx+1:lx+m-1)=y0(lx+1:lx+m-1)-cw(1:m-1); % equivalent to y0=filter(w,1,xx^0); yd=y0(m:end); yd(abs(yd)<eps)=1; nr=length(r); wlx=ones(lx,1); wlxx=ones(lxx,1); y=zeros(lx,nr); mr=max(abs(r)); % max moment to calculate mk=zeros(1,mr); % list of moments mk(-r(r<0))=1; % choose the moments we need to calculate maxr=max(r); if maxr>1 mk(1:maxr)=1; end ml=find(mk>0); lml=length(ml); if lml mlx=mk; mlx(ml)=1:lml; % mapping from moment into ml wlml=ones(1,lml); xm=filter(w,1,xx(:,wlml).^ml(wlxx,:)); % calculate all the moments xm=xm(m:end,:)./yd(:,wlml); % remove the useless start values and normalize end fr=find(r<0); if length(fr) y(:,fr)=xm(:,mlx(-r(fr))); % zero-centred moments end fr=find(r==0); % 0'th moment if length(fr) y(:,fr)=1; end fr=find(r==1); % 1'st moment about mean if length(fr) y(:,fr)=0; end fr=find(r==2); % 2'nd moment about mean if length(fr) yfr=xm(:,2)-xm(:,1).^2; y(:,fr)=yfr(:,ones(1,length(fr))); end if maxr>2 mon=[1 -1]; bc=[1 -2 1]; am=zeros(lx,maxr); am(:,1)=xm(:,1); % copy the mean across ml=2:maxr; wlml=ones(1,length(ml)); am(:,2:end)=xm(:,wlml).^ml(ones(lx,1),:); % calculate powers of the mean for i=3:maxr bc=conv(bc,mon); % calculate binomial coefficients fr=find(r==i); if length(fr) yfr=xm(:,i)+sum(xm(:,i-1:-1:2).*am(:,1:i-2).*bc(wlx,2:i-1),2)+am(:,i)*(bc(i)+bc(i+1)); y(:,fr)=yfr(:,ones(1,length(fr))); end end end